options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # not see parameters str1..., str2,... using regex as exc='[str1,str2,...]' 
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(4,1),mar=c(1,5,1,1))
      drw[,1,] |> plot(type='l',xlab='',ylab=pr)
      drw[,2,] |> plot(type='l',xlab='',ylab=pr)
      drw[,3,] |> plot(type='l',xlab='',ylab=pr)
      drw[,4,] |> plot(type='l',xlab='',ylab=pr)
      par(mar=c(3,5,3,3))
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}

explanatory variable obs.x have noise

xN(x0.sx),yN(b0+b1*x0,s)

ex8-0.stan

normal regression without explanatory variable’s noise

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}

ex9.stan

normal regression with explanatory variable’s noise

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x10;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
  real<lower=0> sx;
  vector[N] x0;
}  
model{
  for(i in 1:N){
    x[i]~normal(x0[i],sx);
    y[i]~normal(b0+b1*x0[i],s);
  }
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x0[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] x1;
  vector[N1] y1;
  for(i in 1:N1){
    x1[i]=normal_rng(x10[i],sx);
    m1[i]=b0+b1*x10[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)

n1=10

#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -9805.7 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       26      -38.3056   0.000609277   0.000161797           1           1       59    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -38.31
##    b0         8.48
##    b1         1.71
##    s          4.12
##    m0[1]     37.20
##    m0[2]     19.52
##    m0[3]     27.93
##    m0[4]      7.43
##    m0[5]     25.25
##    m0[6]     42.66
##    m0[7]     38.88
##    m0[8]     22.16
##    m0[9]     43.71
##    m0[10]    23.40
##    m0[11]     4.86
##    m0[12]    38.17
##    m0[13]    46.06
##    m0[14]    32.24
##    m0[15]    36.19
##    m0[16]    20.36
##    m0[17]    25.05
##    m0[18]    41.02
##    m0[19]    14.27
##    m0[20]    31.03
##    y0[1]     35.78
##    y0[2]     27.80
##    y0[3]     26.60
##    y0[4]      8.80
##    y0[5]     23.02
##    y0[6]     39.06
##    y0[7]     40.91
##    y0[8]     20.69
##    y0[9]     42.88
##    y0[10]    23.26
##    y0[11]     4.09
##    y0[12]    44.86
##    y0[13]    56.49
##    y0[14]    31.82
##    y0[15]    35.52
##    y0[16]    23.46
##    y0[17]    25.41
##    y0[18]    41.26
##    y0[19]    16.71
##    y0[20]    34.24
##    m1[1]      4.86
##    m1[2]      9.44
##    m1[3]     14.02
##    m1[4]     18.59
##    m1[5]     23.17
##    m1[6]     27.75
##    m1[7]     32.33
##    m1[8]     36.91
##    m1[9]     41.48
##    m1[10]    46.06
##    y1[1]      1.90
##    y1[2]      9.00
##    y1[3]     18.24
##    y1[4]     19.46
##    y1[5]     20.95
##    y1[6]     26.51
##    y1[7]     34.63
##    y1[8]     40.97
##    y1[9]     42.02
##    y1[10]    43.58
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -38.58 -38.23 1.38 1.16 -41.38 -37.10 1.01      421      392
##    b0       8.37   8.36 2.29 2.04   4.47  12.03 1.02      272      264
##    b1       1.71   1.72 0.17 0.16   1.45   1.98 1.01      346      426
##    s        4.66   4.56 0.85 0.79   3.54   6.21 1.00     1167     1252
##    m0[1]   37.23  37.24 1.36 1.37  34.97  39.37 1.00     2056     1473
##    m0[2]   19.46  19.48 1.46 1.36  17.07  21.76 1.01      301      283
##    m0[3]   27.91  27.88 1.12 1.09  26.09  29.76 1.01      610      889
##    m0[4]    7.31   7.32 2.38 2.10   3.26  11.15 1.02      272      265
##    m0[5]   25.21  25.19 1.18 1.11  23.26  27.16 1.01      427      528
##    m0[6]   42.70  42.70 1.71 1.63  39.81  45.43 1.00     1270     1119
##    m0[7]   38.91  38.92 1.45 1.47  36.46  41.19 1.00     1811     1410
##    m0[8]   22.11  22.11 1.31 1.22  19.93  24.21 1.01      337      373
##    m0[9]   43.77  43.76 1.79 1.70  40.74  46.66 1.00     1157      965
##    m0[10]  23.35  23.34 1.25 1.18  21.27  25.36 1.01      364      417
##    m0[11]   4.73   4.77 2.61 2.33   0.40   8.92 1.02      273      268
##    m0[12]  38.20  38.21 1.41 1.43  35.82  40.40 1.00     1929     1377
##    m0[13]  46.13  46.11 1.98 1.86  42.78  49.26 1.00      967      923
##    m0[14]  32.24  32.24 1.15 1.14  30.35  34.10 1.00     1420     1308
##    m0[15]  36.21  36.23 1.30 1.32  34.04  38.28 1.00     2123     1361
##    m0[16]  20.30  20.32 1.41 1.32  17.98  22.56 1.01      310      322
##    m0[17]  25.02  25.00 1.18 1.12  23.05  26.97 1.01      418      520
##    m0[18]  41.06  41.06 1.60 1.55  38.38  43.55 1.00     1477     1331
##    m0[19]  14.18  14.19 1.83 1.63  11.08  17.03 1.02      277      221
##    m0[20]  31.02  31.01 1.13 1.12  29.20  32.86 1.00     1107     1253
##    y0[1]   37.23  37.26 4.92 4.60  29.25  45.43 1.00     1907     1810
##    y0[2]   19.33  19.40 4.98 4.86  11.36  27.30 1.00     1668     1759
##    y0[3]   27.79  27.79 4.83 4.66  20.06  35.60 1.00     1786     1643
##    y0[4]    7.35   7.31 5.17 4.85  -0.95  16.11 1.00      784     1357
##    y0[5]   25.28  25.30 4.88 4.74  17.55  33.24 1.00     1765     1709
##    y0[6]   42.57  42.44 4.98 4.69  34.63  50.81 1.00     1884     1679
##    y0[7]   38.74  38.72 5.06 4.78  30.39  46.99 1.00     2097     1950
##    y0[8]   22.12  22.08 4.84 4.59  14.11  30.11 1.00     1622     1879
##    y0[9]   43.74  43.68 5.19 4.85  35.03  52.44 1.00     2040     2010
##    y0[10]  23.51  23.48 4.95 4.51  15.61  31.52 1.00     1702     1872
##    y0[11]   4.71   4.56 5.41 5.26  -4.06  13.67 1.00      975     1533
##    y0[12]  38.24  38.24 4.86 4.64  30.33  46.45 1.00     1881     1917
##    y0[13]  46.05  46.18 5.20 5.03  37.55  54.25 1.00     1897     1798
##    y0[14]  32.21  32.19 4.82 4.48  24.26  40.05 1.00     1936     1929
##    y0[15]  36.23  36.14 4.94 4.87  28.37  44.33 1.00     1789     1870
##    y0[16]  20.39  20.44 4.87 4.77  12.44  28.44 1.00     1782     1929
##    y0[17]  24.98  24.96 5.00 4.82  16.95  33.12 1.00     2138     2035
##    y0[18]  40.88  40.84 4.99 4.82  32.67  48.97 1.00     2110     2010
##    y0[19]  14.22  14.33 4.96 4.68   5.95  22.33 1.00     1364     1541
##    y0[20]  31.03  31.12 4.91 4.96  23.09  38.85 1.00     2118     1994
##    m1[1]    4.73   4.77 2.61 2.33   0.40   8.92 1.02      273      268
##    m1[2]    9.33   9.33 2.21 1.97   5.56  12.82 1.02      272      256
##    m1[3]   13.93  13.94 1.84 1.65  10.79  16.80 1.02      276      223
##    m1[4]   18.53  18.54 1.52 1.41  15.99  20.94 1.01      294      250
##    m1[5]   23.13  23.11 1.26 1.19  21.02  25.15 1.01      358      393
##    m1[6]   27.73  27.70 1.12 1.09  25.90  29.59 1.01      592      848
##    m1[7]   32.33  32.33 1.15 1.15  30.43  34.19 1.00     1444     1303
##    m1[8]   36.93  36.94 1.34 1.36  34.69  39.04 1.00     2092     1520
##    m1[9]   41.53  41.54 1.63 1.56  38.79  44.08 1.00     1411     1298
##    m1[10]  46.13  46.11 1.98 1.86  42.78  49.26 1.00      967      923
##    y1[1]    4.70   4.73 5.39 5.20  -4.18  13.71 1.01      784     1249
##    y1[2]    9.23   9.32 5.12 4.90   0.97  17.57 1.00      955     1429
##    y1[3]   14.03  14.11 5.09 4.80   5.80  22.60 1.00     1164     1455
##    y1[4]   18.49  18.36 5.13 4.90  10.31  26.75 1.00     1459     1828
##    y1[5]   23.15  23.26 4.80 4.61  14.95  30.67 1.00     1788     1743
##    y1[6]   27.84  27.78 4.90 4.64  20.10  36.04 1.00     1454     1597
##    y1[7]   32.30  32.37 4.82 4.61  24.24  40.06 1.00     1833     1971
##    y1[8]   36.86  36.81 4.88 4.77  28.67  44.79 1.00     2054     1972
##    y1[9]   41.45  41.49 5.00 4.91  33.26  49.55 1.00     2108     2012
##    y1[10]  46.16  46.14 5.18 5.00  37.85  54.59 1.00     1872     1976

#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)

mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 2 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 2 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 3 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 4 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 1 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 4 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 3 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 1 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 1 finished in 0.5 seconds.
## Chain 4 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 3 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 3 finished in 0.6 seconds.
## Chain 4 finished in 0.5 seconds.
## Chain 2 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 2 finished in 1.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.7 seconds.
## Total execution time: 1.0 seconds.
seeMCMC(mcmc,'[m,x]')
##  variable   mean median    sd   mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -43.66 -47.24 15.06 11.33 -62.24 -10.23 1.17       17       36
##    b0       8.53   8.48  2.18  2.11   4.99  12.16 1.01      722      962
##    b1       1.70   1.71  0.15  0.15   1.44   1.95 1.01      676     1234
##    s        3.11   3.25  1.63  1.75   0.26   5.59 1.08       52       32
##    sx       1.85   1.92  0.95  1.03   0.33   3.34 1.02       88       94
##    x0[1]   17.02  16.97  1.16  0.94  15.24  18.92 1.01     1702     1370
##    x0[2]    4.09   4.15  2.08  2.52   0.64   6.91 1.02      113      929
##    x0[3]   11.41  11.41  1.19  0.90   9.38  13.30 1.02     2083     1154
##    x0[4]    2.53   2.43  2.50  3.41  -0.89   6.38 1.06       80      369
##    x0[5]    8.43   8.53  1.51  1.68   5.97  10.55 1.01      204     1138
##    x0[6]   19.73  19.78  1.24  1.06  17.66  21.72 1.02     1076     1364
##    x0[7]   18.38  18.25  1.25  1.09  16.49  20.48 1.01      584     1638
##    x0[8]    6.72   6.85  1.52  1.59   4.19   8.85 1.01      241     1054
##    x0[9]   21.48  21.28  1.41  1.23  19.49  23.86 1.01      362     1340
##    x0[10]   7.97   8.09  1.30  1.19   5.76   9.94 1.01      571      982
##    x0[11]  -1.71  -1.79  1.39  1.15  -4.00   0.63 1.01      781     1377
##    x0[12]  18.10  17.94  1.33  1.19  16.15  20.33 1.01      583     1325
##    x0[13]  21.94  21.95  1.30  1.05  19.81  24.08 1.02     1441     1256
##    x0[14]  15.71  15.61  1.67  2.07  13.36  18.37 1.02      119     1485
##    x0[15]  14.80  14.84  1.55  1.79  12.32  17.06 1.01      198     1374
##    x0[16]   6.87   6.92  1.22  0.98   4.78   8.72 1.02     1887     1176
##    x0[17]   9.93   9.93  1.19  0.96   7.96  11.78 1.02     1800     1811
##    x0[18]  20.47  20.31  1.58  1.68  18.28  23.14 1.01      171     1242
##    x0[19]   2.48   2.63  1.45  1.29   0.00   4.55 1.01      556      950
##    x0[20]  12.70  12.79  1.19  1.03  10.64  14.58 1.01     1061     1400
##    m0[1]   37.45  37.55  1.88  1.51  34.26  40.50 1.04     4144     2460
##    m0[2]   15.55  15.52  3.41  4.40  10.90  20.91 1.02      103     1627
##    m0[3]   27.93  27.95  1.96  1.49  24.62  31.13 1.05     3573     1781
##    m0[4]   12.88  13.18  4.37  5.61   5.66  18.63 1.03       87      725
##    m0[5]   22.90  22.85  2.49  2.94  19.38  26.92 1.01      184     1950
##    m0[6]   42.04  41.87  2.08  1.74  38.73  45.49 1.02     1288     2140
##    m0[7]   39.75  39.97  2.06  1.75  36.27  42.98 1.01     1003     2271
##    m0[8]   20.01  19.92  2.49  2.79  16.33  24.02 1.01      212     2618
##    m0[9]   45.00  45.22  2.29  2.03  41.03  48.48 1.01      437     2514
##    m0[10]  22.13  21.98  2.13  2.03  18.81  25.59 1.01      434     2085
##    m0[11]   5.71   5.97  2.41  2.06   1.58   9.49 1.01      944     1799
##    m0[12]  39.28  39.41  2.15  1.90  35.71  42.67 1.01      635     2519
##    m0[13]  45.78  45.68  2.13  1.85  42.20  49.32 1.04     2592     2314
##    m0[14]  35.23  35.33  2.77  3.48  30.80  39.17 1.02      127     2494
##    m0[15]  33.70  33.56  2.65  3.12  29.96  38.01 1.01      158     1579
##    m0[16]  20.25  20.19  1.97  1.63  17.04  23.42 1.05     3162     2536
##    m0[17]  25.44  25.56  1.97  1.57  22.19  28.55 1.03     2369     2122
##    m0[18]  43.28  43.43  2.58  2.92  38.82  47.02 1.01      174     2189
##    m0[19]  12.82  12.62  2.35  2.16   9.23  16.79 1.01      464     1718
##    m0[20]  30.12  29.93  2.00  1.79  26.93  33.41 1.01     1054     2096
##    y0[1]   37.45  37.59  3.99  3.10  30.83  44.02 1.03     3997     3174
##    y0[2]   15.42  14.53  4.82  4.52   9.02  24.25 1.01      211     2523
##    y0[3]   27.82  27.90  4.00  3.02  21.13  34.45 1.03     3706     3046
##    y0[4]   12.94  13.98  5.52  5.76   2.86  19.82 1.02      131     1738
##    y0[5]   22.79  22.05  4.22  3.37  16.89  30.40 1.01      715     2482
##    y0[6]   42.09  41.72  4.15  3.06  35.48  49.21 1.03     2630     1902
##    y0[7]   39.76  40.27  4.07  3.17  32.57  45.93 1.02     3946     3147
##    y0[8]   20.05  19.25  4.32  3.45  13.83  27.89 1.01      843     2497
##    y0[9]   44.98  45.53  4.10  3.26  37.64  51.25 1.01     2032     2933
##    y0[10]  22.13  21.58  4.09  3.23  16.10  29.52 1.02     2703     2862
##    y0[11]   5.75   6.13  4.20  3.23  -1.46  12.58 1.03     2341     1912
##    y0[12]  39.24  39.78  4.16  3.33  31.94  45.63 1.02     2884     2850
##    y0[13]  45.77  45.63  4.09  3.18  38.96  52.49 1.04     4118     3323
##    y0[14]  35.27  36.09  4.52  3.72  27.09  41.53 1.01      416     3069
##    y0[15]  33.81  33.09  4.40  3.65  27.60  41.95 1.01      399     1781
##    y0[16]  20.22  20.15  4.01  3.10  13.46  27.13 1.03     3854     2612
##    y0[17]  25.47  25.74  4.06  3.23  18.56  32.07 1.03     4008     2482
##    y0[18]  43.39  44.16  4.37  3.50  35.31  49.69 1.01      766     3071
##    y0[19]  12.71  12.15  4.24  3.30   6.16  20.24 1.01     2123     2426
##    y0[20]  30.19  29.79  4.02  3.19  23.81  37.15 1.02     3538     3066
##    m1[1]    4.92   4.84  2.46  2.40   0.93   9.05 1.01      703     1036
##    m1[2]    9.48   9.44  2.11  2.03   6.07  13.01 1.01      730     1011
##    m1[3]   14.04  14.00  1.78  1.70  11.20  16.95 1.01      793     1046
##    m1[4]   18.61  18.58  1.48  1.41  16.23  21.05 1.01      922     1164
##    m1[5]   23.17  23.16  1.26  1.19  21.13  25.25 1.01     1162     1199
##    m1[6]   27.74  27.73  1.15  1.10  25.86  29.64 1.00     1441     1335
##    m1[7]   32.30  32.28  1.18  1.17  30.46  34.23 1.00     1499     1231
##    m1[8]   36.87  36.86  1.35  1.31  34.75  39.08 1.00     1233     1138
##    m1[9]   41.43  41.40  1.60  1.58  38.88  44.06 1.00     1023     1183
##    m1[10]  45.99  46.01  1.91  1.89  42.87  49.11 1.00      906     1157
##    x1[1]   -2.14  -2.12  2.07  1.49  -5.65   1.21 1.01     3796     2739
##    x1[2]    0.54   0.55  2.10  1.51  -2.96   4.11 1.01     4025     2613
##    x1[3]    3.24   3.25  2.04  1.48  -0.04   6.56 1.01     3691     2511
##    x1[4]    6.00   5.95  2.10  1.52   2.66   9.49 1.01     3738     2530
##    x1[5]    8.54   8.57  2.07  1.51   5.03  11.94 1.01     4007     3023
##    x1[6]   11.28  11.30  2.04  1.52   7.86  14.71 1.01     3872     3034
##    x1[7]   13.96  13.96  2.08  1.51  10.50  17.44 1.01     4115     2804
##    x1[8]   16.63  16.64  2.05  1.44  13.30  20.09 1.01     3647     2816
##    x1[9]   19.39  19.35  2.12  1.47  15.89  23.04 1.01     4045     2954
##    x1[10]  22.05  22.04  2.09  1.51  18.64  25.52 1.00     3760     2756
##    y1[1]    4.95   4.84  4.27  3.67  -2.07  11.96 1.00     1734     2421
##    y1[2]    9.54   9.46  3.99  3.35   3.15  16.48 1.01     1941     2115
##    y1[3]   13.97  14.00  3.93  3.26   7.47  20.52 1.01     2543     2807
##    y1[4]   18.64  18.64  3.85  3.12  12.24  25.09 1.01     2452     2792
##    y1[5]   23.23  23.19  3.70  2.94  17.21  29.42 1.01     3187     3161
##    y1[6]   27.72  27.76  3.64  2.85  21.70  33.67 1.01     3655     2929
##    y1[7]   32.26  32.27  3.69  2.79  25.99  38.29 1.01     3215     2990
##    y1[8]   36.90  36.85  3.75  3.09  30.65  43.24 1.01     3071     2653
##    y1[9]   41.45  41.40  3.85  3.19  35.25  47.73 1.01     2953     2574
##    y1[10]  46.04  46.02  3.99  3.27  39.44  52.61 1.01     2485     2414

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

outlier

ex10.stan

objective variable have outlier, y~cauchy(b0+b1*x,s)

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~cauchy(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=cauchy_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=cauchy_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)

n1=10

x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -33476 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       26      -32.5138   9.53072e-05   0.000215404      0.8888      0.8888       47    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -32.51
##    b0         6.61
##    b1         1.82
##    s          3.08
##    m0[1]     12.08
##    m0[2]     15.21
##    m0[3]      9.34
##    m0[4]     10.17
##    m0[5]     18.18
##    m0[6]     12.41
##    m0[7]     13.50
##    m0[8]     19.13
##    m0[9]     13.09
##    m0[10]    17.96
##    m0[11]    16.80
##    m0[12]    21.80
##    m0[13]    11.14
##    m0[14]    16.90
##    m0[15]     9.82
##    m0[16]     9.75
##    m0[17]    18.32
##    m0[18]    22.90
##    m0[19]     8.19
##    m0[20]     8.93
##    y0[1]     10.65
##    y0[2]     11.30
##    y0[3]     12.60
##    y0[4]      8.52
##    y0[5]     21.04
##    y0[6]      9.84
##    y0[7]     11.83
##    y0[8]     17.60
##    y0[9]     14.14
##    y0[10]    15.32
##    y0[11]    14.40
##    y0[12]    19.48
##    y0[13]     9.38
##    y0[14]    11.98
##    y0[15]    14.44
##    y0[16]     5.96
##    y0[17]    18.74
##    y0[18]    24.28
##    y0[19]     4.13
##    y0[20]     7.25
##    m1[1]      8.19
##    m1[2]      9.82
##    m1[3]     11.46
##    m1[4]     13.09
##    m1[5]     14.73
##    m1[6]     16.36
##    m1[7]     17.99
##    m1[8]     19.63
##    m1[9]     21.26
##    m1[10]    22.90
##    y1[1]      6.03
##    y1[2]      9.41
##    y1[3]     15.37
##    y1[4]      7.93
##    y1[5]     13.85
##    y1[6]     12.73
##    y1[7]     17.58
##    y1[8]     17.76
##    y1[9]     24.50
##    y1[10]    25.85
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -33.09 -32.73 1.40 1.14 -35.79 -31.58 1.01      439      849
##    b0       6.59   6.63 1.70 1.60   3.77   9.28 1.02      287      348
##    b1       1.82   1.82 0.35 0.34   1.28   2.41 1.02      349      415
##    s        3.52   3.41 0.68 0.60   2.64   4.72 1.00     1493     1305
##    m0[1]   12.06  12.08 0.91 0.90  10.54  13.51 1.01      387      552
##    m0[2]   15.19  15.22 0.81 0.80  13.82  16.51 1.00     1242     1305
##    m0[3]    9.31   9.34 1.26 1.20   7.22  11.26 1.01      299      354
##    m0[4]   10.15  10.17 1.14 1.10   8.26  11.92 1.01      310      403
##    m0[5]   18.16  18.17 1.08 1.04  16.39  19.93 1.00     1259     1095
##    m0[6]   12.39  12.41 0.88 0.88  10.91  13.79 1.01      415      622
##    m0[7]   13.47  13.48 0.81 0.79  12.11  14.81 1.00      579      898
##    m0[8]   19.11  19.11 1.21 1.19  17.11  21.08 1.01     1063     1101
##    m0[9]   13.06  13.08 0.83 0.80  11.66  14.43 1.00      501      749
##    m0[10]  17.94  17.94 1.05 1.02  16.21  19.65 1.00     1311     1267
##    m0[11]  16.78  16.80 0.92 0.90  15.27  18.26 1.00     1547     1321
##    m0[12]  21.78  21.77 1.64 1.61  19.16  24.43 1.01      692      858
##    m0[13]  11.12  11.14 1.01 0.99   9.43  12.70 1.01      337      444
##    m0[14]  16.87  16.89 0.93 0.90  15.34  18.37 1.00     1536     1321
##    m0[15]   9.80   9.82 1.19 1.13   7.80  11.64 1.01      304      382
##    m0[16]   9.73   9.75 1.20 1.15   7.72  11.58 1.01      303      382
##    m0[17]  18.30  18.32 1.10 1.06  16.50  20.10 1.00     1228     1105
##    m0[18]  22.87  22.86 1.83 1.79  19.99  25.83 1.01      624      868
##    m0[19]   8.17   8.21 1.44 1.37   5.80  10.42 1.02      289      318
##    m0[20]   8.91   8.94 1.32 1.28   6.71  10.98 1.01      295      339
##    y0[1]   11.94  11.92 3.69 3.49   5.87  18.10 1.00     1801     1883
##    y0[2]   15.14  15.15 3.69 3.50   9.06  20.98 1.00     2109     1819
##    y0[3]    9.30   9.48 3.77 3.67   3.24  15.37 1.00     1700     1930
##    y0[4]   10.15  10.02 3.75 3.56   4.29  16.53 1.00     1741     2060
##    y0[5]   18.13  18.18 3.76 3.55  11.85  24.28 1.00     1907     1750
##    y0[6]   12.40  12.43 3.64 3.38   6.35  18.22 1.00     1672     1734
##    y0[7]   13.37  13.36 3.74 3.47   7.20  19.72 1.00     1875     1869
##    y0[8]   19.14  19.13 3.86 3.79  12.82  25.23 1.00     1816     1907
##    y0[9]   13.09  13.11 3.75 3.66   6.93  19.23 1.00     2012     1839
##    y0[10]  18.07  18.16 3.72 3.51  11.94  23.95 1.00     1872     1767
##    y0[11]  16.82  16.79 3.67 3.48  10.86  22.77 1.00     1822     2023
##    y0[12]  21.82  21.69 4.03 3.75  15.40  28.51 1.00     1650     1545
##    y0[13]  10.96  10.93 3.73 3.66   4.92  16.95 1.00     1729     1800
##    y0[14]  16.88  16.78 3.73 3.58  10.68  23.09 1.00     2014     1683
##    y0[15]   9.62   9.69 3.87 3.53   3.49  16.02 1.00     1255     1667
##    y0[16]   9.60   9.65 3.85 3.76   3.32  15.80 1.00     1349     1515
##    y0[17]  18.30  18.28 3.80 3.68  12.19  24.34 1.00     1808     1931
##    y0[18]  23.08  23.15 4.06 3.81  16.39  29.64 1.00     1381     1884
##    y0[19]   8.13   8.19 3.88 3.60   1.60  14.59 1.00     1282     1608
##    y0[20]   8.99   8.92 3.88 3.77   2.77  15.33 1.00     1264     1872
##    m1[1]    8.17   8.21 1.44 1.37   5.80  10.42 1.02      289      318
##    m1[2]    9.80   9.83 1.19 1.13   7.80  11.64 1.01      304      382
##    m1[3]   11.44  11.46 0.98 0.96   9.80  12.97 1.01      350      469
##    m1[4]   13.07  13.08 0.83 0.80  11.67  14.44 1.00      502      749
##    m1[5]   14.71  14.73 0.80 0.77  13.38  16.00 1.00     1038     1114
##    m1[6]   16.34  16.36 0.88 0.86  14.84  17.75 1.00     1553     1343
##    m1[7]   17.97  17.98 1.06 1.03  16.24  19.69 1.00     1302     1220
##    m1[8]   19.61  19.59 1.29 1.24  17.51  21.70 1.01      969     1011
##    m1[9]   21.24  21.23 1.55 1.53  18.75  23.76 1.01      739      898
##    m1[10]  22.87  22.86 1.83 1.79  19.99  25.83 1.01      624      868
##    y1[1]    8.13   8.18 3.77 3.52   1.98  14.48 1.00     1209     1768
##    y1[2]    9.80   9.82 3.81 3.57   3.35  15.99 1.00     1404     1615
##    y1[3]   11.48  11.45 3.68 3.43   5.50  17.72 1.00     1631     1833
##    y1[4]   13.04  13.00 3.68 3.54   7.03  19.04 1.00     1570     1796
##    y1[5]   14.72  14.65 3.71 3.52   8.60  20.87 1.00     1921     1515
##    y1[6]   16.34  16.35 3.67 3.39  10.37  22.43 1.00     1959     1892
##    y1[7]   17.99  18.12 3.76 3.46  11.77  24.09 1.00     1930     1812
##    y1[8]   19.59  19.54 3.82 3.72  13.38  25.77 1.00     2054     2105
##    y1[9]   21.43  21.43 3.89 3.68  15.07  27.86 1.00     1658     1780
##    y1[10]  22.89  22.95 4.07 3.91  16.11  29.36 1.00     1536     1588

mdl=cmdstan_model('./ex10.stan') 
fn(mdl,data)
## Initial log joint probability = -82.5218 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       20      -6.62665   0.000293134   0.000533178           1           1       37    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__      -6.63
##    b0         5.52
##    b1         1.89
##    s          0.37
##    m0[1]     11.20
##    m0[2]     14.45
##    m0[3]      8.35
##    m0[4]      9.22
##    m0[5]     17.53
##    m0[6]     11.53
##    m0[7]     12.66
##    m0[8]     18.51
##    m0[9]     12.24
##    m0[10]    17.29
##    m0[11]    16.09
##    m0[12]    21.28
##    m0[13]    10.22
##    m0[14]    16.19
##    m0[15]     8.85
##    m0[16]     8.78
##    m0[17]    17.67
##    m0[18]    22.42
##    m0[19]     7.16
##    m0[20]     7.93
##    y0[1]     11.16
##    y0[2]     14.58
##    y0[3]      6.20
##    y0[4]     10.80
##    y0[5]     18.49
##    y0[6]     13.01
##    y0[7]     12.37
##    y0[8]     17.91
##    y0[9]     12.47
##    y0[10]    17.90
##    y0[11]    16.49
##    y0[12]    21.19
##    y0[13]    10.17
##    y0[14]    16.21
##    y0[15]     8.86
##    y0[16]     8.47
##    y0[17]    17.78
##    y0[18]    21.65
##    y0[19]     7.03
##    y0[20]    23.14
##    m1[1]      7.16
##    m1[2]      8.86
##    m1[3]     10.55
##    m1[4]     12.25
##    m1[5]     13.94
##    m1[6]     15.64
##    m1[7]     17.33
##    m1[8]     19.03
##    m1[9]     20.72
##    m1[10]    22.42
##    y1[1]      6.38
##    y1[2]     15.73
##    y1[3]      2.45
##    y1[4]     11.74
##    y1[5]     15.35
##    y1[6]     14.99
##    y1[7]     17.36
##    y1[8]     19.29
##    y1[9]     20.33
##    y1[10]    22.61
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable  mean median     sd  mad     q5   q95 rhat ess_bulk ess_tail
##    lp__   -9.08  -8.74   1.27 1.00 -11.52 -7.73 1.01      458      618
##    b0      5.52   5.53   0.29 0.28   5.05  5.98 1.00      696      616
##    b1      1.90   1.90   0.07 0.06   1.80  2.01 1.00      821      859
##    s       0.45   0.43   0.14 0.13   0.26  0.72 1.01      749      603
##    m0[1]  11.21  11.21   0.15 0.15  10.98 11.48 1.00      927     1059
##    m0[2]  14.47  14.46   0.16 0.15  14.23 14.75 1.00     1928     1421
##    m0[3]   8.36   8.36   0.21 0.21   8.03  8.70 1.00      722      745
##    m0[4]   9.23   9.22   0.19 0.19   8.93  9.53 1.00      750      799
##    m0[5]  17.56  17.54   0.22 0.20  17.23 17.95 1.00     1697     1536
##    m0[6]  11.55  11.55   0.15 0.14  11.32 11.81 1.00      978     1050
##    m0[7]  12.68  12.67   0.15 0.14  12.46 12.94 1.00     1278     1365
##    m0[8]  18.55  18.53   0.25 0.23  18.18 18.99 1.00     1550     1377
##    m0[9]  12.26  12.25   0.15 0.14  12.03 12.51 1.00     1135     1255
##    m0[10] 17.32  17.31   0.21 0.20  17.01 17.70 1.00     1737     1488
##    m0[11] 16.12  16.10   0.19 0.17  15.84 16.45 1.00     1879     1560
##    m0[12] 21.32  21.31   0.33 0.31  20.84 21.91 1.00     1282     1252
##    m0[13] 10.23  10.22   0.17 0.16   9.96 10.51 1.00      815      861
##    m0[14] 16.22  16.20   0.19 0.18  15.94 16.55 1.00     1869     1542
##    m0[15]  8.86   8.86   0.20 0.20   8.55  9.18 1.00      736      784
##    m0[16]  8.79   8.79   0.20 0.20   8.47  9.11 1.00      734      782
##    m0[17] 17.70  17.69   0.22 0.21  17.37 18.10 1.00     1673     1511
##    m0[18] 22.46  22.44   0.36 0.35  21.92 23.10 1.00     1218     1175
##    m0[19]  7.17   7.17   0.24 0.24   6.78  7.56 1.00      702      671
##    m0[20]  7.94   7.94   0.22 0.22   7.59  8.30 1.00      712      713
##    y0[1]  11.38  11.19  23.69 0.65   8.16 14.31 1.00     1883     1828
##    y0[2]  14.39  14.45  12.30 0.67  11.65 17.67 1.00     2075     1933
##    y0[3]   3.49   8.36 232.58 0.69   5.62 11.42 1.00     1700     1685
##    y0[4]   9.28   9.22  12.85 0.71   6.00 12.45 1.00     1900     1728
##    y0[5]  17.92  17.58   9.71 0.73  14.95 20.25 1.00     1763     1735
##    y0[6]  12.26  11.54  45.29 0.66   7.72 15.07 1.00     2050     1849
##    y0[7]  15.00  12.67 146.49 0.69   9.46 15.20 1.00     2092     1833
##    y0[8]  18.65  18.54   8.23 0.72  15.60 21.82 1.00     1815     1599
##    y0[9]  11.56  12.25  49.36 0.69   9.40 15.56 1.00     1859     1972
##    y0[10] 16.92  17.32  52.59 0.68  14.19 20.02 1.00     1920     1950
##    y0[11] 15.83  16.10   8.68 0.66  13.05 18.79 1.00     1779     2020
##    y0[12] 21.82  21.32   8.55 0.79  18.86 24.49 1.00     1916     1881
##    y0[13] 12.31  10.19 115.81 0.70   7.24 13.12 1.00     2002     2041
##    y0[14] 16.09  16.21  18.35 0.72  13.01 18.89 1.00     2026     1920
##    y0[15]  9.37   8.86  11.58 0.70   5.92 11.60 1.00     1695     1860
##    y0[16]  6.87   8.78  72.88 0.70   6.13 11.48 1.00     2018     1857
##    y0[17] 18.18  17.71  30.30 0.68  14.76 20.70 1.00     1879     1913
##    y0[18] 20.34  22.42 102.03 0.78  19.74 25.42 1.00     2145     1955
##    y0[19]  7.08   7.20   7.43 0.69   4.48  9.54 1.00     2062     1943
##    y0[20]  8.16   7.94  15.75 0.73   4.77 11.07 1.00     1770     1907
##    m1[1]   7.17   7.17   0.24 0.24   6.78  7.56 1.00      702      671
##    m1[2]   8.87   8.86   0.20 0.20   8.55  9.18 1.00      736      784
##    m1[3]  10.56  10.56   0.16 0.16  10.31 10.84 1.00      844      855
##    m1[4]  12.26  12.26   0.15 0.14  12.04 12.52 1.00     1135     1255
##    m1[5]  13.96  13.95   0.15 0.15  13.73 14.22 1.00     1790     1399
##    m1[6]  15.66  15.65   0.18 0.17  15.39 15.98 1.00     1931     1455
##    m1[7]  17.36  17.34   0.22 0.20  17.04 17.74 1.00     1730     1488
##    m1[8]  19.06  19.04   0.26 0.24  18.67 19.53 1.00     1489     1365
##    m1[9]  20.76  20.75   0.31 0.29  20.30 21.32 1.00     1321     1304
##    m1[10] 22.46  22.44   0.36 0.35  21.92 23.10 1.00     1218     1175
##    y1[1]   7.25   7.16  23.55 0.69   4.13 10.02 1.00     1618     1669
##    y1[2]   8.38   8.83  11.89 0.68   5.50 11.46 1.00     1966     1955
##    y1[3]  10.46  10.57   9.25 0.68   7.47 13.45 1.00     1992     1981
##    y1[4]  11.99  12.25  33.92 0.67   9.65 14.76 1.00     2023     1853
##    y1[5]  14.60  13.95  26.01 0.72  11.26 17.28 1.00     2011     1932
##    y1[6]  15.01  15.64  19.49 0.74  12.14 18.33 1.00     1690     1891
##    y1[7]  17.08  17.35  24.06 0.67  14.43 20.13 1.00     1824     1710
##    y1[8]  19.39  19.03  23.02 0.71  16.06 21.74 1.00     1872     1964
##    y1[9]  21.07  20.74  16.73 0.76  18.10 24.02 1.00     1901     1991
##    y1[10] 22.34  22.43  12.73 0.79  18.90 25.53 1.00     1791     1745

censored data

objective variable has NA

(x,y) i=1-N
(x0,y0) i=1-N0
x1 i=1-N1, y1=NA
(x,y)~N((mx,my),(sx2,sy2,sxy))
(x0,y0)~N((mx,my),(sx2,sy2,sxy))
x1~N(mx,sx2)

ex11-0.stan

data{
  int N0;
  array[N0] vector[2] xy;
  int N1;
  vector[N1] x1;
}
parameters{
  vector[2] m;
  cov_matrix[2] s;
}
model{
  target+=multi_normal_lpdf(xy | m, s);
  x1~normal(m[1],s[1,1]^.5);
}
generated quantities{
  vector[2] xy1;
  xy1=multi_normal_rng(m,s);
  real cr;
  cr=s[1,2]/(s[1,1]*s[2,2])^.5;
}
n=30
x=runif(n,0,9)
y=rnorm(n,10+3*x,4)
cor(x,y)
## [1] 0.769
qplot(x,y)

L=4
n0=sum(x>L)
x0=x[x>L]
y0=y[x>L]
x1=x[x<=L]
n1=sum(x<=L)
cor(x0,y0)
## [1] 0.597
qplot(x0,y0)

mdl=cmdstan_model('./ex11-0.stan') 

data=list(N0=n0,xy=matrix(c(x0,y0),ncol=2),N1=n1,x1=x1)

mle=mdl$optimize(data=data)
## Initial log joint probability = -154798 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Exception: multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite.  last conditional variance is 0. (in '/tmp/Rtmp3zdloC/model-bd3037f1494a.stan', line 12, column 2 to column 39) 
## Exception: multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite.  last conditional variance is 0. (in '/tmp/Rtmp3zdloC/model-bd3037f1494a.stan', line 12, column 2 to column 39) 
##       55      -135.468   0.000281555   0.000319153           1           1       98    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__    -135.47
##    m[1]       6.06
##    m[2]      27.70
##    s[1,1]     4.35
##    s[2,1]    10.70
##    s[1,2]    10.70
##    s[2,2]    50.79
##    xy1[1]     6.59
##    xy1[2]    34.60
##    cr         0.72
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.3 seconds.
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 0.3 seconds.
## Chain 3 finished in 0.3 seconds.
## Chain 4 finished in 0.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.3 seconds.
## Total execution time: 0.4 seconds.
seeMCMC(mcmc,ch=F)
##  variable    mean  median    sd   mad      q5     q95 rhat ess_bulk ess_tail
##    lp__   -131.26 -130.93  1.71  1.58 -134.62 -129.19 1.00      597      934
##    m[1]      6.05    6.04  0.43  0.44    5.35    6.76 1.00      886      938
##    m[2]     27.71   27.78  1.66  1.64   24.84   30.26 1.02      457      586
##    s[1,1]    5.74    5.39  1.82  1.52    3.53    9.11 1.00     1362     1548
##    s[2,1]   14.04   12.84  6.54  5.41    6.02   26.41 1.00      464      631
##    s[1,2]   14.04   12.84  6.54  5.41    6.02   26.41 1.00      464      631
##    s[2,2]   70.26   63.60 29.40 23.50   36.55  127.07 1.00      502      735
##    xy1[1]    6.13    6.12  2.48  2.40    2.11   10.16 1.00     1762     1703
##    xy1[2]   27.77   27.80  8.84  8.38   12.90   41.95 1.00     1710     1708
##    cr        0.69    0.72  0.14  0.13    0.42    0.87 1.01      460      797

xy=mcmc$draws('xy1')
cor(xy[,,1],xy[,,2])
## [1] 0.707
qplot(xy[,,1],xy[,,2])

objective variable is censored

y i=1-N, y~N(m,s)  
  actual          ya i=1-Na
  lower censored  yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
  upper censored  yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)

cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))

P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)

ex11-1.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  real L;
  int Nl;
  vector[Nl] xl;
  real U;
  int Nu;
  vector[Nu] xu;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
  for(i in 1:Nl)
    target+=normal_lcdf(L | b0+b1*xl[i],s);
  for(i in 1:Nu)
    target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)

xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
mdl=cmdstan_model('./ex8-0.stan')

data=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)

mle=mdl$optimize(data=data)
## Initial log joint probability = -3560.18 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       20      -23.3087    0.00437841   0.000468356           1           1       37    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -23.31
##    b0        17.68
##    b1         1.48
##    s          3.64
##    m0[1]     24.66
##    m0[2]     21.36
##    m0[3]     23.03
##    m0[4]     26.06
##    m0[5]     25.95
##    m0[6]     20.02
##    m0[7]     21.26
##    m0[8]     22.64
##    m0[9]     21.12
##    m0[10]    25.26
##    m0[11]    20.54
##    m0[12]    24.27
##    m0[13]    28.13
##    y0[1]     18.42
##    y0[2]     13.38
##    y0[3]     20.01
##    y0[4]     26.45
##    y0[5]     25.72
##    y0[6]     19.16
##    y0[7]     25.18
##    y0[8]     27.03
##    y0[9]     25.98
##    y0[10]    22.96
##    y0[11]    18.84
##    y0[12]    24.86
##    y0[13]    29.96
##    m1[1]     17.86
##    m1[2]     19.25
##    m1[3]     20.65
##    m1[4]     22.04
##    m1[5]     23.44
##    m1[6]     24.84
##    m1[7]     26.23
##    m1[8]     27.63
##    m1[9]     29.02
##    m1[10]    30.42
##    y1[1]     13.67
##    y1[2]     21.14
##    y1[3]     21.89
##    y1[4]     28.13
##    y1[5]     28.91
##    y1[6]     29.89
##    y1[7]     25.31
##    y1[8]     28.42
##    y1[9]     23.86
##    y1[10]    30.50
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m')
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -23.64 -23.29 1.31 1.07 -26.18 -22.18 1.00      585      809
##    b0      17.51  17.35 3.15 3.09  12.76  22.97 1.01      309      285
##    b1       1.53   1.56 0.74 0.74   0.25   2.68 1.00      361      494
##    s        4.46   4.29 1.10 0.96   3.06   6.38 1.00     1402     1344
##    m0[1]   24.72  24.69 1.39 1.32  22.52  27.03 1.00     1797     1326
##    m0[2]   21.32  21.24 1.62 1.51  18.79  24.02 1.01      378      516
##    m0[3]   23.03  23.02 1.26 1.14  20.97  25.06 1.00      780      929
##    m0[4]   26.17  26.19 1.82 1.67  23.27  29.10 1.00     1097     1151
##    m0[5]   26.05  26.07 1.78 1.64  23.19  28.92 1.00     1144     1137
##    m0[6]   19.93  19.81 2.12 2.00  16.71  23.44 1.01      324      333
##    m0[7]   21.21  21.14 1.66 1.54  18.64  23.98 1.01      370      486
##    m0[8]   22.63  22.59 1.31 1.18  20.49  24.76 1.01      605      786
##    m0[9]   21.06  20.97 1.70 1.59  18.41  23.94 1.01      361      491
##    m0[10]  25.34  25.36 1.55 1.44  22.86  27.90 1.00     1540     1199
##    m0[11]  20.46  20.34 1.91 1.78  17.52  23.68 1.01      337      364
##    m0[12]  24.32  24.32 1.32 1.23  22.20  26.49 1.00     1822     1367
##    m0[13]  28.30  28.33 2.67 2.57  23.84  32.56 1.00      672      931
##    y0[1]   24.75  24.76 4.82 4.64  16.87  32.84 1.00     2079     1973
##    y0[2]   21.16  21.12 4.85 4.52  13.48  28.91 1.00     1539     1578
##    y0[3]   22.75  22.81 4.68 4.41  15.17  30.32 1.00     1961     1910
##    y0[4]   26.18  26.28 4.99 4.75  17.96  34.30 1.00     1690     1822
##    y0[5]   26.10  26.11 4.79 4.48  18.73  33.92 1.00     1966     1822
##    y0[6]   20.06  19.93 5.01 4.83  12.13  28.07 1.00      864     1370
##    y0[7]   21.10  21.10 4.90 4.57  13.37  28.95 1.00     1248     1766
##    y0[8]   22.48  22.37 4.66 4.37  14.91  29.90 1.00     1967     1971
##    y0[9]   21.03  21.12 4.89 4.52  13.16  28.96 1.00     1030     1753
##    y0[10]  25.33  25.15 4.87 4.47  17.59  33.66 1.00     1979     1846
##    y0[11]  20.38  20.21 4.94 4.59  12.54  28.50 1.00     1060     1492
##    y0[12]  24.21  24.17 4.99 4.61  16.30  32.43 1.00     1953     1884
##    y0[13]  28.28  28.50 5.33 4.91  19.29  36.87 1.00     1343     1538
##    m1[1]   17.69  17.55 3.07 3.01  13.07  23.05 1.01      309      278
##    m1[2]   19.13  18.99 2.44 2.37  15.41  23.18 1.01      315      290
##    m1[3]   20.58  20.46 1.87 1.75  17.68  23.74 1.01      340      375
##    m1[4]   22.02  21.96 1.43 1.28  19.73  24.35 1.01      459      681
##    m1[5]   23.46  23.45 1.25 1.13  21.40  25.55 1.00     1102     1270
##    m1[6]   24.90  24.90 1.43 1.36  22.63  27.28 1.00     1766     1326
##    m1[7]   26.34  26.38 1.88 1.74  23.34  29.34 1.00     1033     1165
##    m1[8]   27.78  27.81 2.45 2.34  23.73  31.71 1.00      727      958
##    m1[9]   29.23  29.27 3.08 2.91  24.06  34.00 1.00      605      868
##    m1[10]  30.67  30.69 3.73 3.56  24.37  36.53 1.00      531      740
##    y1[1]   17.66  17.56 5.53 5.12   8.79  26.98 1.01      826     1236
##    y1[2]   19.18  19.11 5.32 5.02  10.69  27.95 1.00      960     1510
##    y1[3]   20.73  20.72 5.03 4.51  12.68  29.07 1.00      999     1256
##    y1[4]   22.06  22.04 4.99 4.49  14.23  30.00 1.00     1480     1974
##    y1[5]   23.46  23.47 4.61 4.39  16.30  30.75 1.00     2163     2037
##    y1[6]   25.01  24.92 4.78 4.37  17.23  32.81 1.00     2107     1934
##    y1[7]   26.27  26.42 5.07 4.79  18.12  34.59 1.00     1601     1777
##    y1[8]   27.68  27.61 5.32 5.05  19.18  36.34 1.00     1430     1856
##    y1[9]   29.21  29.28 5.50 5.13  20.39  37.90 1.00     1213     1407
##    y1[10]  30.70  30.77 5.95 5.59  20.97  40.36 1.00      879     1461

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

data=list(N=n,x=xya$x,y=xya$y,
          L=L,Nl=nl,xl=xyl$x,
          U=U,Nu=nu,xu=xyu$x,
          N1=n1,x1=x1)
mdl=cmdstan_model('./ex11-1.stan') 

mle=mdl$optimize(data=data)
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Initial log joint probability = -117.14 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       31      -29.0816   0.000876035   4.91276e-07           1           1       41    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -29.08
##    b0        12.65
##    b1         2.78
##    s          4.30
##    m0[1]     25.78
##    m0[2]     19.58
##    m0[3]     22.71
##    m0[4]     28.41
##    m0[5]     28.20
##    m0[6]     17.06
##    m0[7]     19.39
##    m0[8]     21.98
##    m0[9]     19.12
##    m0[10]    26.91
##    m0[11]    18.03
##    m0[12]    25.05
##    m0[13]    32.30
##    y0[1]     25.50
##    y0[2]     26.64
##    y0[3]     25.35
##    y0[4]     31.36
##    y0[5]     24.37
##    y0[6]     20.58
##    y0[7]     15.76
##    y0[8]     21.35
##    y0[9]     16.90
##    y0[10]    29.00
##    y0[11]    17.16
##    y0[12]    30.25
##    y0[13]    34.52
##    m1[1]     12.99
##    m1[2]     15.61
##    m1[3]     18.24
##    m1[4]     20.86
##    m1[5]     23.48
##    m1[6]     26.11
##    m1[7]     28.73
##    m1[8]     31.35
##    m1[9]     33.98
##    m1[10]    36.60
##    y1[1]     12.40
##    y1[2]      9.25
##    y1[3]     23.55
##    y1[4]     22.33
##    y1[5]     30.85
##    y1[6]     33.42
##    y1[7]     19.85
##    y1[8]     37.26
##    y1[9]     25.63
##    y1[10]    39.91
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,'m',ch=T)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -29.25 -28.90 1.36 1.08 -31.80 -27.79 1.01      564      616
##    b0      11.71  11.91 3.13 3.01   6.40  16.34 1.02      241      322
##    b1       3.04   2.98 0.71 0.68   2.02   4.25 1.01      277      437
##    s        5.33   5.04 1.40 1.15   3.63   8.03 1.00     1099     1073
##    m0[1]   26.07  25.98 1.41 1.32  23.91  28.51 1.00     1554     1214
##    m0[2]   19.29  19.38 1.71 1.56  16.51  21.95 1.01      319      402
##    m0[3]   22.71  22.69 1.35 1.25  20.43  24.90 1.01      586      803
##    m0[4]   28.94  28.82 1.76 1.63  26.28  32.06 1.00     1350     1029
##    m0[5]   28.72  28.60 1.72 1.59  26.11  31.76 1.00     1394     1120
##    m0[6]   16.53  16.63 2.17 2.04  12.93  19.81 1.02      262      351
##    m0[7]   19.08  19.16 1.74 1.60  16.22  21.78 1.01      312      400
##    m0[8]   21.91  21.92 1.40 1.27  19.58  24.17 1.01      477      717
##    m0[9]   18.78  18.87 1.78 1.62  15.86  21.51 1.01      303      381
##    m0[10]  27.30  27.20 1.53 1.42  24.99  30.02 1.00     1624     1233
##    m0[11]  17.59  17.70 1.98 1.87  14.31  20.59 1.01      276      364
##    m0[12]  25.27  25.20 1.36 1.26  23.18  27.47 1.00     1339     1262
##    m0[13]  33.20  33.02 2.52 2.31  29.52  37.93 1.00      624      970
##    y0[1]   25.86  25.78 5.84 5.42  16.32  35.41 1.00     2094     1931
##    y0[2]   19.17  19.26 5.73 5.32   9.81  28.21 1.00     1464     1673
##    y0[3]   22.59  22.61 5.68 5.42  13.44  31.64 1.00     1894     1856
##    y0[4]   29.18  28.96 5.96 5.51  19.57  38.99 1.00     1753     1962
##    y0[5]   28.78  28.86 5.85 5.57  19.90  38.24 1.00     1818     1921
##    y0[6]   16.63  16.73 6.12 5.43   6.71  26.64 1.00     1357     1450
##    y0[7]   19.33  19.24 5.87 5.22  10.02  28.80 1.01     1261     1528
##    y0[8]   22.02  21.96 5.70 5.24  12.90  31.12 1.00     1746     1664
##    y0[9]   18.67  18.59 5.92 5.42   9.67  28.19 1.00     1559     1601
##    y0[10]  27.37  27.24 5.56 5.19  18.51  36.55 1.00     1758     1860
##    y0[11]  17.58  17.75 5.95 5.26   7.49  26.91 1.00     1459     1560
##    y0[12]  25.17  25.29 5.55 5.10  16.27  34.38 1.00     2022     2014
##    y0[13]  33.40  33.36 6.08 5.90  23.70  43.19 1.00     1287     1678
##    m1[1]   12.08  12.28 3.05 2.93   6.89  16.59 1.02      241      321
##    m1[2]   14.95  15.10 2.47 2.36  10.84  18.61 1.02      250      319
##    m1[3]   17.82  17.91 1.94 1.83  14.59  20.77 1.01      281      364
##    m1[4]   20.69  20.74 1.52 1.37  18.23  23.12 1.01      379      518
##    m1[5]   23.56  23.51 1.32 1.22  21.41  25.69 1.00      760     1087
##    m1[6]   26.43  26.34 1.44 1.35  24.23  28.95 1.00     1617     1235
##    m1[7]   29.30  29.20 1.81 1.68  26.56  32.50 1.00     1288     1039
##    m1[8]   32.17  32.02 2.32 2.12  28.78  36.48 1.00      726     1054
##    m1[9]   35.04  34.82 2.89 2.69  30.83  40.38 1.00      521      940
##    m1[10]  37.90  37.62 3.50 3.29  32.78  44.27 1.00      440      791
##    y1[1]   12.01  12.25 6.23 5.83   1.40  21.80 1.00      707     1196
##    y1[2]   14.91  15.39 6.06 5.60   4.34  23.96 1.00      989     1122
##    y1[3]   17.84  17.85 5.69 5.27   8.50  27.11 1.00     1458     1842
##    y1[4]   20.69  20.70 5.56 5.09  11.59  29.64 1.00     1884     1878
##    y1[5]   23.53  23.59 5.66 5.14  14.23  32.81 1.00     1870     1958
##    y1[6]   26.42  26.40 5.71 5.20  17.35  35.91 1.00     1954     1926
##    y1[7]   29.15  29.15 5.93 5.39  19.40  39.06 1.00     1952     1957
##    y1[8]   32.13  31.89 5.89 5.40  22.96  41.91 1.00     1656     1848
##    y1[9]   34.95  34.72 6.25 5.79  25.34  45.69 1.00     1652     1561
##    y1[10]  37.84  37.51 6.70 6.29  27.67  48.79 1.00     1169     1487

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

sensitivity/specificity analysis

ex14.stan

estimating sens and spec

data {
  int N;
  array[N] int x;
  array[N] int y;
}
parameters {
  real<lower=0,upper=1> p;
  real<lower=0,upper=1> se;
  real<lower=0,upper=1> sp;
}
model {
  p~uniform(0,1);
  se~uniform(0,1);
  sp~uniform(0,1);
  for (i in 1:N) {
    y[i]~bernoulli(x[i]*se+(1-x[i])*(1-sp));
  }
}
generated quantities {
  real ppv;
  real npv;
  ppv=se*p/((se*p)+((1-p)*(1-sp)));
  npv=(1-p)*sp/(((1-p)*sp)+(p*(1-se)));
}
n=20
data=list(N=n,
          x=sample(0:1,n,replace=T),
          y=sample(0:1,n,replace=T))

mdl=cmdstan_model('./ex14.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -13.0211 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5      -12.0405   0.000152385   1.84742e-05      0.9347      0.9347        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -12.04
##      p        0.70
##      se       0.75
##      sp       0.37
##      ppv      0.73
##      npv      0.40
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -18.13 -17.81 1.32 1.12 -20.67 -16.70 1.00      708      896
##      p      0.49   0.48 0.29 0.37   0.04   0.95 1.00      673      642
##      se     0.71   0.73 0.12 0.12   0.50   0.89 1.00     2566     1340
##      sp     0.41   0.40 0.14 0.15   0.18   0.65 1.00     2449     1411
##      ppv    0.52   0.53 0.29 0.37   0.05   0.95 1.01      685      667
##      npv    0.57   0.62 0.29 0.36   0.07   0.97 1.00      734      579

2x2 cross table

Effect occur y=1 with probabilty p01, p11 from each cause x{0,1}
event frequncy nxy of effect y{0,1} by cause x{0,1}
n01~B(n0.,p0)
n11~B(n1.,p1)

n01=n0p0, n00=n0(1-p0)
n11=n1p1, n10=n1(1-p1)

p00=n00/n=n0(1-p0)/n, p01=n01/n=n0p0/n
p10=n10/n=n1(1-p1)/n, p11=n11/n=n1p1/n

Cramer'V  (chi2/n/(min(row,column)-1))^.5
  in 2x2  
  crv =(n11*n00-n10*n01)/(n0.*n1.*n.0*n.1)^.5
      =(n0(1-p0)n1p1-n0p0n1(1-p1))/(n0n1(n0(1-p0)+n1(1-p1))(n0p0n1p1))^.5

kappa coefficient   k=(po-pe)/(1-pe)
  po: Observed agreement (proportion of times both raters agreed)
  pe: Expected agreement under independence
      po=p00+p11
        =(n0(1-p0)+n1p1)/n
      pe=(p00+p01)(p00+p10)(p10+p11)(p01+p11)
        =n0/n*(n0(1-p0)+n1(1-p1))/n*(n0p0+n1p1)/n*n1/n

ex16-1.stan

data {
  int n;
  int n0;
  int n01;
  int n1;
  int n11;
}
parameters {
  real<lower=0,upper=1> p0;
  real<lower=0,upper=1> p1;
}
model {
  n01~binomial(n0,p0);
  n11~binomial(n1,p1);
}
generated quantities {
  real RR;
  RR=p1/p0;
  real OR;
  OR=(p1/(1-p1))/(p0/(1-p0));
}

ex16-2.stan

data {
  int n;
  int n0;
  int n01;
  int n1;
  int n11;
}
parameters {
  real<lower=0,upper=1> p0;
  real<lower=0,upper=1> p1;
}
model {
  n01~binomial(n0,p0);
  n11~binomial(n1,p1);
}
generated quantities {
  real RR;
  RR=p1/p0;
  real OR;
  OR=(p1/(1-p1))/(p0/(1-p0));
  real crv;
  crv=(n0*(1-p0)*n1*p1-n0*p0*n1*(1-p1))/(n0*n1*(n0*(1-p0)+n1*(1-p1))*(n0*p0+n1*p1))^.5;
  real k;
  real po;
  po=(n0*(1-p0)+n1*p1)/n;
  real pe;
  pe=n0/n*(n0*(1-p0)+n1*(1-p1))/n*(n0*p0+n1*p1)/n*n1/n;
  k=(po-pe)/(1-pe);
}
n0=30
n01=rbinom(1,n0,0.3)
n1=30
n11=rbinom(1,n1,0.6)
data=list(n=n0+n1,n0=n0,n01=n01,n1=n1,n11=n11)

mdl=cmdstan_model('./ex16-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -71.8484 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5      -39.2858    0.00286964   0.000356069      0.9698      0.9698        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -39.29
##      p0       0.33
##      p1       0.60
##      RR       1.80
##      OR       3.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -43.24 -42.93 1.04 0.72 -45.32 -42.27 1.00      942     1027
##      p0     0.34   0.34 0.08 0.09   0.21   0.48 1.00     1966     1351
##      p1     0.59   0.59 0.08 0.09   0.46   0.73 1.00     2150     1342
##      RR     1.85   1.75 0.60 0.52   1.11   2.91 1.00     1986     1453
##      OR     3.34   2.88 1.98 1.55   1.23   6.95 1.00     1976     1473

data=list(n=n0+n1,n0=n0,n01=n01,n1=n1,n11=n11)

mdl=cmdstan_model('./ex16-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -47.5796 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        4      -39.2858   0.000782113   1.21406e-05           1           1        7    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -39.29
##      p0       0.33
##      p1       0.60
##      RR       1.80
##      OR       3.00
##      crv      0.27
##      k        0.61
##      po       0.63
##      pe       0.06
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -43.24 -42.93 1.04 0.72 -45.32 -42.27 1.00      942     1027
##      p0     0.34   0.34 0.08 0.09   0.21   0.48 1.00     1966     1351
##      p1     0.59   0.59 0.08 0.09   0.46   0.73 1.00     2150     1342
##      RR     1.85   1.75 0.60 0.52   1.11   2.91 1.00     1986     1453
##      OR     3.34   2.88 1.98 1.55   1.23   6.95 1.00     1976     1473
##      crv    0.25   0.26 0.12 0.13   0.05   0.44 1.00     1976     1486
##      k      0.60   0.60 0.06 0.07   0.49   0.70 1.00     1975     1473
##      po     0.62   0.63 0.06 0.06   0.53   0.72 1.00     1976     1450
##      pe     0.06   0.06 0.00 0.00   0.06   0.06 1.00     1519     1277

bimodal distribution, mixed normal distribution

n=100
y0=rnorm(n,0,1)
y1=rnorm(n,-5,1)
y2=rnorm(n,5,1)
y=sample(c(y0,y1,y2),n)
density(y) |> plot()

EM algorithm

library(mclust)

rst=Mclust(y)
summary(rst)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust E (univariate, equal variance) model with 4 components: 
## 
##  log-likelihood   n df  BIC  ICL
##            -230 100  8 -497 -500
## 
## Clustering table:
##  1  2  3  4 
## 36 24  7 33
rst$parameters
## $pro
## [1] 0.3600 0.2415 0.0664 0.3321
## 
## $mean
##     1     2     3     4 
## -5.32 -0.14  3.09  5.44 
## 
## $variance
## $variance$modelName
## [1] "E"
## 
## $variance$d
## [1] 1
## 
## $variance$G
## [1] 4
## 
## $variance$sigmasq
## [1] 0.51
plot(rst)

ex17-1.stan

data {
  int N;
  vector[N] y;;
}
parameters {
  simplex[3] h; //ratio of mix
  real m1;
  real m2;
  real m3;
  real<lower=0> s1;
  real<lower=0> s2;
  real<lower=0> s3;
}
model {
  s1~cauchy(0,5);
  s2~cauchy(0,5);
  s3~cauchy(0,5);
  for (i in 1:N) {
    vector[3] p;
    p[1]=log(h[1]) + normal_lpdf(y[i] | m1, s1);
    p[2]=log(h[2]) + normal_lpdf(y[i] | m2, s2);
    p[3]=log(h[3]) + normal_lpdf(y[i] | m3, s3);
    target+=log_sum_exp(p);
  }
}
mdl=cmdstan_model('./ex17-1.stan')

data=list(N=n,y=y)
mle=mdl$optimize(data=data)
## Initial log joint probability = -474.65 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       93      -289.352   0.000914839    0.00664515      0.8877      0.8877      137    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__  -289.35
##      h[1]     0.92
##      h[2]     0.05
##      h[3]     0.03
##      m1       0.08
##      m2      -0.38
##      m3       0.31
##      s1       4.77
##      s2       0.08
##      s3       0.03
mcmc=goMCMC(mdl,data,smp=2000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:    1 / 2100 [  0%]  (Warmup) 
## Chain 1 Iteration:  101 / 2100 [  4%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:    1 / 2100 [  0%]  (Warmup) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:    1 / 2100 [  0%]  (Warmup) 
## Chain 3 Iteration:  101 / 2100 [  4%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:    1 / 2100 [  0%]  (Warmup) 
## Chain 1 Iteration: 1100 / 2100 [ 52%]  (Sampling) 
## Chain 3 Iteration: 1100 / 2100 [ 52%]  (Sampling) 
## Chain 2 Iteration:  101 / 2100 [  4%]  (Sampling) 
## Chain 1 Iteration: 2100 / 2100 [100%]  (Sampling) 
## Chain 4 Iteration:  101 / 2100 [  4%]  (Sampling) 
## Chain 1 finished in 1.3 seconds.
## Chain 3 Iteration: 2100 / 2100 [100%]  (Sampling) 
## Chain 3 finished in 1.4 seconds.
## Chain 4 Iteration: 1100 / 2100 [ 52%]  (Sampling) 
## Chain 2 Iteration: 1100 / 2100 [ 52%]  (Sampling) 
## Chain 4 Iteration: 2100 / 2100 [100%]  (Sampling) 
## Chain 4 finished in 28.3 seconds.
## Chain 2 Iteration: 2100 / 2100 [100%]  (Sampling) 
## Chain 2 finished in 31.6 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 15.7 seconds.
## Total execution time: 31.8 seconds.
seeMCMC(mcmc,ch=F)
##  variable    mean  median     sd   mad      q5     q95 rhat ess_bulk ess_tail
##      lp__ -257.32 -260.91  15.38 23.12 -276.31 -239.76 2.02        5       35
##      h[1]    0.13    0.10   0.13  0.14    0.00    0.32 1.73        6      128
##      h[2]    0.35    0.35   0.06  0.06    0.26    0.44 1.14       18       56
##      h[3]    0.51    0.52   0.15  0.21    0.31    0.73 1.81        5       43
##      m1    614.14    1.77 863.80  4.51   -0.30 2522.67 2.17        5       13
##      m2     -0.04   -0.27   5.29  7.76   -5.50    5.50 1.90        5       36
##      m3      0.10    0.50   4.17  6.10   -5.43    5.30 2.84        4       26
##      s1     16.09    1.10 388.25  0.71    0.58   30.95 1.41       10      159
##      s2      0.84    0.81   0.21  0.17    0.51    1.22 1.43        8       28
##      s3      2.06    1.93   1.18  1.59    0.70    3.83 2.16        5       30